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Abstract 

In compressible turbulent combustion/nonequilibrium flows, the constructions of numeri- 
cal schemes for (a) stable and accurate simulation of turbulence with strong shocks, and 
(b) obtaining correct propagation speed of discontinuities for stiff reacting terms on “coarse 
grids” share one important ingredient - minimization of numerical dissipation while main- 
taining numerical stability. Here “coarse grids” means standard mesh density requirement 
for accurate simulation of typical non-reacting flows. This dual requirement to achieve both 
numerical stability and accuracy with zero or minimal use of numerical dissipation is most 
often conflicting for existing schemes that were designed for non-reacting flows. The goal 
of this paper is to relate numerical dissipations that are inherited in a selected set of high or- 
der shock-capturing schemes with the onset of wrong propagation speed of discontinuities 
for two representative stiff detonation wave problems. 

Key words: High order numerical methods, Numerical methods for turbulence with 
shocks, Stiff source terms, Wrong propagation speed of discontinuities 


1. Introduction 

To make accurate predictions in compressible (magnetized) turbulent combus- 
tion/nonequilibrium flows, one has to deal with the equations that describe time-dependent 
non-equilibrium effects, combustion, advanced thermodynamic models, and magnetic 
fields (MHD). Numerical simulation is challenging because of the conflicting requirements 
for numerical methods to be accurate enough to resolve the small scales of turbulence but 
robust enough to handle shock waves without generating spurious numerical noise. In ad- 
dition, the different physics models have different time scales that, when underresolved, 
might interact numerically to produce erroneous results. Furthermore, the appearance of 
the source terms in modeling flow problems containing finite-rate chemistry or combus- 
tion poses additional numerical difficulties beyond that for solving non-reacting turbulent 
flows. The so-called stiff source term problem [7] is a well-known example. For stiff reac- 
tions it is well known that the wrong propagation speed of discontinuities occurs due to the 
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under-resolved numerical solutions in both space and time. Schemes to improve the wrong 
propagation speed of discontinuities for systems of stiff reacting flows remain a challenge 
for algorithm development [15]. 

In addition to the minimization of numerical dissipation while maintaining numerical 
stability in compressible turbulence with strong shocks, Yee & Sjogreen and Yee & Sweby 
[19, 20, 16] discussed a general framework for the design of such schemes. Yee & Sjogreen 
[22], Wang et al. [15] and references cited therein present their recent progress on the 
subject. In [25] a short overview of this recent progress is given. Two very important 
numerical challenges are “Stiffness and Nonlinearity of Source Terms”. 

The objective of the present paper is to a gain a deeper understanding on the behavior 
of four high order shock-capturing schemes with the onset of wrong propagation speed of 
discontinuities for two representative stiff detonation wave problems. The test cases consist 
of the Arrhenius ID Chapman- Jouguet (C-J) detonation wave [4, 13] and a 2D Heaviside 
detonation wave [1], These are the same two test cases considered in [15]. The considered 
four schemes are the fifth-order WENO, “WEN05”, the newly developed subcell resolu- 
tion version of WEN05, “WEN05/SR” [15], the Yee & Sjogreen nonlinear filter version 
of WEN05 using a local flow sensor to further limit the amount of WEN05 numerical 
dissipation, “WEN05fi”, and the Durcros split version of WEN05fi, “WEN05fi+split” 
[21, 22], All of the four methods use the Roe’s average states. For the temporal discretiza- 
tion the classical fourth-order Runge-Kutta method (RK4) is used. See the aforementioned 
references for the development of these schemes. 

WEN05/SR [15] is a newly developed modified fractional step method which solves the 
convection step and reaction step separately. In the convection step any high order shock- 
capturing method can be used. In the reaction step an ODE solver is applied, but with the 
computed flow variables in the shock region modified by the Harten subcell resolution idea 

[3]. 

WEN05fi is the filter version of WEN05. On the first stage a full time step by RK4 
is performed. For this stage the sixth-order central spatial base scheme is used. On the 
second stage the solution is filtered by the dissipative portion of WEN05 in conjunction 
with a local wavelet flow sensor [22], The wavelet flow sensor indicates the locations and 
the amount where shock-capturing dissipations are needed and leaves the remaining region 
free of numerical dissipation contamination. WEN05fi+split is WEN05fi applied to the 
Ducros et al. split form of the governing equation [2] before the application of WEN05fi. 
The Ducros et al. split form is a preprocessing step to condition the governing equation(s) 
before the application of high order central schemes. This preprocessing step improves 
numerical stability and is widely used in numerical modeling and simulation of turbulent 
flows. 

The comparison of the performance of the four schemes is largely based on the degree 
that each method captures the correct location and jump size of the stiff reaction front for 
coarse grids. Here “coarse grids” means standard mesh density requirement for accurate 
simulation of typical non-reacting flows. It is remarked that, in order to resolve the sharp 
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reaction zone, sufficiently many grid points in this zone are still needed. The behavior of 
these schemes in the vicinity of a sharp reaction zone with several levels of grid refinement 
will be briefly touched upon. 

2. Numerical methods 

Consider a 2D reactive Euler equation with two chemical states of burnt gas and unbumt 
gas and a single irreversible reaction. Without heat conduction and viscosity, the system can 
be written as 


Pt + (pu) x + ( pv)y 

= 0 

( 1 ) 

(, pu)t + ( pU 2 + //),. + (pUV)y 

= 0 

( 2 ) 

(, PV)t + ( PUV) X + ( pV 2 + p)y 

= 0 

( 3 ) 

Et + (u(E + p)) x + (v(E + p))y 

= 0 

( 4 ) 

(, pz)t + (, puz) x + ( pvz) y 

= —K(T)pz, 

( 5 ) 


where p(x, y, t ) is the mixture density, u(x, y. t ) and v(x, y, t ) are the mixture x- and y- 
velocities, E(x, y, t ) is the mixture total energy per unit volume, p(x, y, t ) is the pressure, 
z(x, y, t ) is the mass fraction of the unbumt gas, K(T) is the chemical reaction rate and 
T(x, y, t) is the temperature. The pressure is given by 

P={ 7 - 1)(£ - ^p(u 2 + v 2 ) - q 0 pz ), (6) 

where the temperature T =~ p and q 0 is the chemical heat released in the reaction. 

The reaction rate K(T) is modeled by an Arrhenius law 

K(T) = K 0 exp(^^j, (7) 

where K 0 is the reaction rate constant and T ign is the ignition temperature. The reaction 
rate may be also modeled in the Heaviside form 

K(T ) = { 1//£ T ~ Tian (8) 

[ ’ 1 0 T<T ign ’ W 

where e is the reaction time and 1/e is roughly equal to K 0 . 

Here only the newly developed high order finite difference method with subcell resolu- 
tion for advection equations with stiff source terms [15] in 2D is briefly summarized. The 
key aspect of the filter counterpart of the WENO schemes are included at the end of the 
section. 
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2.1. High order finite difference methods with subcell resolution for advection equations 
with stiff source terms 

The general fractional step approach based on Strang-splitting [12] for equation 

U t + F(U) X + G(U) y = S(U ) (9) 

is as follows. The numerical solution at time level t n+ 1 is approximated by 


V n + 1 = W R{At)A ( — )U\ 


The reaction operator R is over a time step At and the convection operator A is over 
At/2. The two half-step reaction operations over adjacent time steps can be combined to 
save cost. The convection operator A is defined to approximate the solution of the homo- 
geneous part of the problem on the time interval, i.e., 

Ut + F{U) X + G(U) y = 0, t n < t < /„- i . (11) 

The reaction operator R is defined to approximate the solution on a time step of the reaction 
problem: 

^ = S(U), t n < t < t n+ 1 . (12) 

Here, the convection operator consists of, e.g.., WEN05 with Roe flux and RK4 for 
time discretization. If there is no smearing of discontinuities in the convection step, any 
ODE solver can be used as the reaction operator. However, all the standard shock-capturing 
schemes will produce a few transition points in the shock when solving the convection 
equation. These transition points are usually responsible for causing incorrect numerical 
results in the stiff case. Thus we cannot directly apply a standard ODE solver at these 
transition points. Here the Harten’s subcell resolution technique in the reaction step is 
employed. The general idea is as follows. If a point is considered a transition point of 
the shock, information from its neighboring points which are deemed not transition points 
will be used instead. In 2D case we apply the subcell resolution procedure dimension by 
dimension. Here, U T = ip, pu, pv, E, pz). The algorithm proceeds as follows. 

(1) Use a “shock indicator” to identify cells in which discontinuities are believed to be 
situated. One can use any indicator suitable for the particular problem. Here the minmod- 
based shock indicator in [3, 11] is considered. Identify troubled cell I t j in both x- and 
//-directions by applying the shock indicator to, e.g., the mass fraction z. Define the cell I tJ 
as troubled in the ^-direction if | sf \ > |.sf_ 1;/ | and | sf \ > |.sf +1 j| with at least one strict 
inequality, where 

Sij = minmod{:ij + i j - z t] ,z t] - z^j}. (13) 

Similarly we can define the cell I tJ as troubled in the //-direction if | s/ 3 \ > \s/ J _ 1 \ and 
141 > 1 4 + i I with at least one strict inequality where 

4 = minmodj> !J + 1 - z t] ,z t] - z hJ _ ] }. (14) 
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If Iij is only troubled in one direction, we apply the subcell resolution along this direc- 
tion. If Iij is troubled in both directions, we choose the direction which has a larger jump. 
Namely, if |s? | > | s y tJ | , subcell resolution is applied along the ^-direction, otherwise it is 
done along the ^-direction. In the following steps (2)-(3), without loss of generality, we 
assume the subcell resolution is applied in the ^-direction. Assuming I l3 is troubled in the 
^-direction, we apply subcell resolution along the ^-direction. 

In a troubled cell identified above, we continue to identify its neighboring cells. For 
example, we can define Ii+ij as troubled if |sf +1 -| > |,s^L lj and |,sf +lj | > |s? +2 -| and 
similarly define h-ij as troubled if sf_ , ■ > |sf_ 2j | and |sf_ x -| > |sf +1 -|. If the cell 
Ii- S ,j and the cell h+ r ,j (s, r > 0) are the first good cells from the left and the right (i.e., 
Ii-s+i,j and h+r-ij are still troubled cells), we compute the fifth-order ENO interpolation 
polynomials pi- s j(x ) and p i+r j(x) for the cells U- S ,j and Ii +r j, respectively. 

(2) Modify the point values z, tJ , Tij and p i:) in the troubled cell I tJ by the ENO interpo- 
lation polynomials 


f Z ij Pi—s,jiAii Z )i 

\ ~*i] Pi+r,j (*U> z) i 


Tij Pi— s,j i T') , pij Pi—gj^XifP'), if 6 ^ Xi 

Tij Pi+r,] (-A 1 T) , Pij pi+r,j ('A ^ p) ; if 0 < .1', 


(15) 


where the location 9 is determined by the conservation of energy E 



Pi- S ,j(x; E)dx + 



Pi +r j(x ; E)dx = EijAx. 


(16) 


Under certain conditions, it can be shown that there is a unique 9 satisfying Eq. (16), which 
can be solved using, for example, a Newton’s method. If there is no solution for 9 or there is 
more than one solution, we choose z l} = z i+r j, T\ 3 = T l+r] and p VJ = p l+rj . For particular 
problems one can choose any other suitable method for the reconstruction. 

(3) Use Uij instead of U l:1 in the ODE solver if the cell I tJ is a troubled cell. For 
simplicity, explicit Euler is used as the ODE solver. 


(P^)p +1 = (P z )ij + A tS(Tij, Pij , Zij )• 


(17) 


In general, a regular CFL=0.1 can be used in the proposed scheme to produce a stable 
solution. But the solution is very coarse in the reaction zone because of the underresolved 
mesh in time. In order to obtain more accurate results in the reaction zone, we evolve one 
reaction step via N r sub steps, i.e.. 


u n+1 = A 





( 18 ) 


in some numerical examples studied in [15]. 
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2.2. Well-Balanced High Order Filter Schemes for Reacting Flows [9, 21, 22, 24 ] 

Before the application of a high order non-dissipative spatial base scheme, the pre- 
processing step to improve stability had split inviscid flux derivatives of the governing 
equation(s) in the following three ways, depending on the flow types and the desire for 
rigorous mathematical analysis or physical argument. 

• Entropy splitting of Olsson & Oliger [8] and Yee et al. [18, 19]: The resulting form is 
non-conservative and the derivation is based on entropy norm stability with numerical 
boundary closure for the initial value boundary problem. 

• The system form of the Ducros et al. splitting [2]: This is a conservative splitting and 
the derivation is based on physical arguments. 

• Tadmor entropy conservation formulation for systems (Sjogreen & Yee [10]): The 
derivation is based on mathematical analysis. It is a generalization of Tadmor’s en- 
tropy formulation to systems and has not been fully tested on complex flows. 

After the application of a non-dissipative high order spatial base scheme on the split 
form of the governing equation(s), to further improve nonlinear stability from the non- 
dissipative spatial base scheme, the post-processing step of Yee & Sjogreen [21, 22], 
Sjogreen & Yee [9] nonlinearly filtered the solution by a dissipative portion of a high order 
shock-capturing scheme with a local flow sensor. The flow sensor provides locations and 
amounts of built-in shock-capturing dissipation that can be further reduced or eliminated. 
The idea of these nonlinear filter schemes for turbulence with shocks is that, instead of 
solely relying on very high order high-resolution shock-capturing methods for accuracy, the 
filter schemes [17, 18, 9, 21] take advantage of the effectiveness of the nonlinear dissipation 
contained in good shock-capturing schemes as stabilizing mechanisms (a post-processing 
step) at locations where needed. The nonlinear dissipative portion of a high-resolution 
shock-capturing scheme can be any shock-capturing scheme. By design, the flow sensors, 
spatial base schemes and nonlinear dissipation models are standalone modules. Unlike stan- 
dard shock-capturing and/or hybrid shock-capturing methods, the nonlinear filter method 
requires one Riemann solve per dimension per time step, independent of time discretiza- 
tions. The nonlinear filter method is more efficient than its shock-capturing method coun- 
terparts employing the same order of the respective methods. See [22, 24] for the recent 
improvements of the work [17, 18, 9, 21] that are suitable for a wide range of flow speed 
with minimal tuning of scheme parameters. For all the computations shown, the Ducros 
et al. splitting is employed. This is due to the fact that for the subject test cases we need 
a robust conservative splitting as the preprocessing step. The subcell resolution approach 
using the fractional step procedure can carry over to the aforementioned filter schemes as 
well. Some attributes of the high order filter approach are: 

• Spatial Base Scheme: High order and conservative (no flux limiter or Riemann 
solver) 
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Figure 1: Pressure and density (mass fraction of unbumt gas) comparison among four high order shock- 
capturing methods for the C-J detonation problem, Arrhenius case at t = 1.8 using 50 uniform grid points. 


• Physical Viscosity: Contribution of physical viscosity, if it exists, is automatically 
taken into consideration by the base scheme in order to minimize the amount of nu- 
merical dissipation to be used by the filter step 

• Efficiency: One Riemann solve per dimension per time step, independent of time dis- 
cretizations (less CPU time and fewer grid points than their standard shock-capturing 
scheme counterparts) 

• Accuracy: Containment of numerical dissipation via a local wavelet flow sensor 

• Well-balanced scheme: These nonlinear filter schemes are well-balanced schemes for 
certain chemical reacting flows [14] 

• Stiff Combustion with Discontinuities: For some stiff reacting flow test cases the high 
order filter scheme is able to obtain the correct propagation speed of discontinuities, 
whereas the standard high order shock-capturing (e.g., WENO) schemes cannot (see 
the result below) 

• Parallel Algorithm: Suitable for most current supercomputer architectures 


3. Numerical Examples 

The behavior of the considered four methods is investigated on two test cases that were 
considered in [15], The test cases consist of the Arrhenius ID C-J detonation wave and a 2D 
Heaviside detonation wave. Note that the computed solutions by WEN05 and WEN05/SR 
presented here could be slightly different from the results presented in [15] due to the minor 
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Figure 2: Pressure comparison among four high order shock-capturing methods for the C-J detonation prob- 
lem, Arrhenius case at t = 1.8 using 800 uniform grid points. 

differences in the formulation of the problem (governing equation; e.g., different choice of 
variables) and the use of cell centered and cell- vertex formulation of numerical schemes. 

3.1. Chapman- Jouguet ( C-J) Detonation Wave ( Arrhenius Case ) 

The test case is the ID C-J detonation wave (Arrhenius case) [4, 13]. The initial values 
consist of totally burnt gas on the left-hand side and totally unburnt gas on the right-hand 
side. The density, velocity, and pressure of the unburnt gas are given by p u = 1, u u = 0 
and p u = 1. 

The initial state of the burnt gas is calculated from C-J condition: 

Pt=-!>+(!> 2 -c) 1/2 , 

i>uW(' r + i) - p»] 

Pb = , 

7 Pb 

Scj = [ PuU u + ( 7 PbPbY /2 }/ Pu, 

u b = S C j - ( iPb/rhob ) 11 2 , 

where 

b = -p u ~ puQoil ~ 1), (23) 

c = P« + 2 (7- 1Kmo/(7 + 1). (24) 

The heat release q 0 = 25 and the ratio of specific heats is set to 7 = 1.4. The ig- 
nition temperature T ign = 25 and K 0 = 164180. The computation domain is [0,30]. 
Initially, the discontinuity is located at x = 10. At time t = 1.8, the detonation wave 
has moved to x = 22.8. The reference solution is computed by the regular WEN05 


(19) 

( 20 ) 

( 21 ) 

( 22 ) 
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scheme with 10000 uniform grid points and CFL=0.05. Figure 1 shows the pressure and 
mass fraction comparison among the standard WEN05 scheme, WEN05/SR, WEN05fi 
and WEN05fi+split using 50 uniform grid points. For this particular problem and grid 
size, WEN05fi+split compares well with WEN05/SR for the computed pressure solution. 
WEN05/SR and WEN05fi+split can capture the correct structure using fewer grid points 
than those in Helzel et al. and Tosatto &Vigevano [4, 13]. A careful examination of the 50 
coarse grid mass fraction solutions indicates that WEN05fi+split is one grid point ahead 
of WEN05/SR at the discontinuity location when compared to the reference solution. The 
reference solution is obtained by WEN05 using 10, 000 grid points. Since WEN05fi+split 
is less dissipative than WEN05, the restriction of the shock-capturing dissipation using 
the wavelet flow sensor helps to improve the wrong propagation speed of discontinuities 
without the subcell resolution procedure. Figure 2 shows a grid refinement in the hope of 
resolving the narrow reaction zone using 800 uniform grid points. It is interesting to see 
that all of the methods (except WEN05) produce oscillatory solutions in the vicinity of the 
reaction front. This behavior prompted us to perform a systematic six levels of uniform grid 
refinements (200, 400, 800, 1600, 3200 and 6400). As the number of grid point increases, 
this oscillatory behavior in the vicinity of the reaction front becomes more pronounced. It 
appears that for this particular test case the new subcell resolution scheme WEN05/SR, 
WENOfi and WEN05fi+split only perform well with coarse grids. However, for the more 
dissipative scheme WEN05, as we refine the grid, the computed solution gets closer and 
closer to the reference solution. 


3.2. 2D detonation waves 

This example is taken from [1], The chemical reaction is modeled by the Heaviside 
form with the parameters 


7 = 1.4, q 0 = 0.5196 x 10 10 , 


- = 0.5825 x 10 10 , 

£ 


T ign = 0.1155 x 10 10 


in CGS units. Consider a two-dimensional channel of width 0.005 with solid walls at the 


upper and lower boundaries. The computational domain is [0, 0.025] x [0, 0.005]. The initial 


conditions are 


(p,u,v,p, z) 


(p b ,u b ,Q,p b ,Q), if x<£(y), 
(p u ,u u ,0,p u ,l), if x>£(y), 


(25) 


where 


£(y) 


0.004 \y - 0.00251 > 0.001, 

0.005 - | y- 0.0025| \y - 0.0025| < 0.001, 


(26) 


and u v = 0, p u = 1.201 x 10 -3 , p u = 8.321 x 10 5 and u b = 8.162 x 10 4 . Values of p b and 
Ph are defined by Eq. (19) and (20). In this case u b is greater than defined by Eq. (22). 

One important feature of this solution is the appearance of triple points, which travel 
in the transverse direction and reflect from the upper and lower walls. A discussion of 
the mechanisms driving this solution is given in [6], Figures 3 and 4 show the density 
comparison among the standard WEN05 scheme, WEN05/SR and WEN05fi+split using 
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Figure 3: Density computed for the 2D detonation problem at £ = 0.3x 10“ 7 by different methods. From left 
to right: reference solution by the standard WEN05 method using 2000 x 400 uniform grid points, WEN05, 
WEN05/SR and WEN05fi+split using 500 x 100 uniform grid points. 
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Figure 4: Density computed for the 2D detonation problem at£ = 1.7xl0 _7 by different methods. From left 
to right: reference solution by the standard WEN05 method using 2000 x 400 uniform grid points, WEN05, 
WEN05/SR and WEN05fi+split using 500 x 100 uniform grid points. 
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Figure 5: ID cross-section of density at t = 1.7x10 7 by four high order shock-capturing methods for the 2D 
detonation problem using 200 x 40 uniform grid points with the left figure in the vicinity of the discontinuity. 


500 x 100 uniform grid points at two different times. Figure 5 shows the density comparison 
among the standard WEN05 scheme, WEN05/SR, WEN05fi and WEN05fi+split using 
200 x 40 and 500 x 100 uniform grid points. The reference solutions are computed by 
standard WEN05 with 2000 x 400 grid points. Again, WEN05/SR and WEN05fi+split 
are able to obtain the correct shock speed with similar accuracy. WEN05fi gives a slight 
oscillatory solution near x = 0.004. WEN05 and WEN05/SR produce no oscillations at 
the same location. Further improvement of the flow sensor of the filter scheme is needed in 
order to remove the spurious oscillations. Furthermore, for the 500 x 100 grid, WEN05fi 
also obtained the correct shock speed. 

In summary, we demonstrated that the filter version of the WEN05 in conjunction 
with the Ducros et al. splitting (WEN05fi+split) is able to obtain the correct propaga- 
tion speed of discontinuities for two detonation problems. From the result WEN05/SR and 
WEN05fi+split are able to obtain the correct shock speed with similar accuracy, whereas 
this is not the case for WEN05 WEN05fi using the same coarse grids. Using its original 
form [22] without further modification, the accuracy of WEN05fi+split is nearly as good 
as the proposed high-order finite difference schemes with subcell resolution. The next step 
is to examine the subcell resolution version of WEN05fi and WEN05fi+split, and their 
seventh-order counterparts.. 
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